Simple excision of a black hole in 3+1 numerical relativity 
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We describe a simple implementation of black hole excision in 3+1 numerical relativity. We ap- 
ply this technique to a Schwarzschild black hole with octant symmetry in Eddington-Finkelstein 
coordinates and show how one can obtain accurate, long-term stable numerical evolutions. 
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The simulation of a black hole inspiral collision is one of 
the most important open problems facing numerical rel- 
ativity. Traditional techniques using singularity avoiding 
slicings will not be able to follow such a collision since 
problems associated with the stretching of the slice typ- 
ically cause simulations to crash or to become extremely 
inaccurate in time scales far shorter than the orbital time 
scale. Black hole excision techniques (also known as "ap- 
parent horizon boundary condition" appear to be 
the most promising way of eliminating the problem of the 
slice stretching, thus in principle allowing numerical sim- 
ulations to follow the inspiral from well separated holes 
through the merger and the ring-down phase. 

Black hole excision was first attempted successfully by 
Seidel and Suen in spherical symmetry [ 1 1 , and was later 
studied in more detail by Anninos et.al. 2]. However, the 
original idea is older, and Thornburg |^,| has attributed 
it to a suggestion by Unruh from 1984. The idea con- 
sists of two parts: First, one places a boundary inside 
the black hole and excises its interior from the computa- 
tional domain. Second, one uses a shift vector that keeps 
the horizon roughly in the same coordinate location dur- 
ing the evolution ( "horizon tracking" , see [Q ) . Since no 
information can leave the interior of the black hole, exci- 
sion should have no effect on the physics outside. Ideally, 
one would like to know the position of the event horizon 
which marks the true causal boundary, but the global 
character of its definition means that in principle one can 
only locate it once the whole evolution of the spacetime 
is known. The apparent horizon, on the other hand, can 
be located on every time slice and is guaranteed to be in- 
side the event horizon. In practice one therefore needs to 
find the apparent horizon and excise a region contained 
inside it. 

Though black hole excision has been successful in 
spherical symmetry [p]Pj5|-p^ , it has been difficult to 
implement with a 3+1 approach in three-dimensions 
(3D) |]ll|-|r^, where instabilities typically plague the evo- 
lutions (but some progress has been made, sec [ |l5| , p^ ). 
Black hole excision using a characteristic formulation, on 
the other hand, has been very successful in 3D, allowing 
stable evolutions of perturbed black holes for thousands 
of Af's [|l^. However, such characteristic formulations 
are likely to have problems with the development of caus- 
tics in the case of extremely distorted or colliding black 



holes, so the search for a stable 3+1 excision implemen- 
tation is still of great importance. 

Here we present a 3+1 approach to black hole excision 
in 3D that has allowed us to obtain long-term stable, ac- 
curate evolutions of a single black hole spacetime. These 
results are currently limited to simulations in octant sym- 
metry as discussed below. 

I. SIMPLE BLACK HOLE EXCISION 

Though conceptually simple, black hole excision in 3D 
is a complicated problem numerically. First, one has to 
cut a hole in the computational domain that has a spher- 
ical topology and is therefore not well adapted to the 
Cartesian coordinates typically used. Second, one has 
to apply some condition at the boundary of the excised 
region that is stable and respects the causality of the 
physical system. As the excised region is inside a black 
hole, no boundary condition should be needed since all 
the information required to update the boundary comes 
from outside the excised region. However, achieving this 
"boundary without a boundary condition" (BWBC) [00] 
in 3D is difficult, particularly if one uses a formulation of 
the evolution equations that is not hyperbolic. The way 
this problem is usually approached is by using "causal 
differencing" or "causal reconnection" |T^, where 

the computational molecules are adapted to follow the 
causal structure. The mixture of these issues makes it 
difficult in practice to identify what particular element 
of an algorithm is responsible for causing a numerical 
simulation to go unstable. 

In our approach we have simplified the algorithm as 
much as possible, separating out what we believe is es- 
sential to the excision problem. Our algorithm is based 
on the following simplifications: 

• Excise a region adapted to Cartesian coordinates, 
i.e. excise a cube contained inside the horizon. 

• Do not attempt to fulfill the BWBC ideal, and use 
instead a simple but stable boundary condition at 
the excision boundary. 

• Do not use causal differencing. Use instead cen- 
tered differences in all terms except the advection 
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terms on the shift (terms that look like ). For 
these terms use upwind along the shift direction (we 
use the standard ID second-order upwind stencil in 
each of the Cartesian coordinate directions based 
on the sign of the corresponding shift component 
at each point). This is very important, as it is the 
only place where any information about causality 
(i.e. the direction of the shift) enters our scheme. 
Using a centered approximation for these terms re- 
sults in an unstable scheme. 

One can worry that excising a cube will introduce ar- 
tifacts into the evolution, but as long as the boundary 
condition used at the sides of the cube is consistent those 
artifacts will converge away with increased resolution. 
Similarly, one can argue that applying a boundary con- 
dition instead of using causal differencing is inconsistent 
with the physics, but since this condition is applied well 
inside the horizon, any error introduced is unlikely to 
propagate outside the hole. 

II. STATIC BLACK HOLE SPACETIME 

As the first test of our excision algorithm we have 
considered a single static black hole written in "3+1 
Eddington-Finkelstein" (EF) coordinates. These 3-1-1 EF 
coordinates are a simple transformation of the standard 
ingoing EF coordinates |19| to a 3+1 form. The result- 
ing metric has no coordinate singularities, penetrates the 
event horizon, reaches the physical singularity, and is 
manifestly time independent. This makes it ideal for ex- 
cision tests where one can excise the physical singularity 
and try to keep the numerical evolution stable and close 
to static. The 3+1 EF metric has the form 
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jijtTK/3), and = ^^''T'jf. (note that 



ds^ = - (1 - 2M/r) dt^ + (4Af/r) dtdr 
+ {l + 2M/r) dr^+r^dfl^ . 



(1) 



with M the black hole mass and the solid angle el- 
ement. From this metric one can read the values of the 
3-metric, lapse and shift. The extrinsic curvature can 
then be obtained in a straightforward way. 

III. EVOLUTION EQUATIONS 

Formulation — We comment briefly on the formulation 
used for the simulations described below. Our simula- 
tions have been performed using a formulation of the 
3+1 evolution equations developed by Baumgarte and 
Shapiro |Q (^^L based on previous work of Shibata 
and Nakamura |^l[ (SN). The motivation for using this 
BSSN formulation comes from the fact that it has shown 
remarkable stability properties when compared to the 
Arnowitt-Deser-Misner (ADM) formulation ||2^ in a wide 
range of numerical simulations (2^j2^ ^ . 

The BSSN variables are defined in terms of the 
spatial metric 7^ and the extrinsic curvature Kij 
as: ^ = ln(det7y)/12, 7y=e"'''^7y, trK = ^''^Kij, 



for the explicit form of 
28 1 for an analysis that in- 



det.g = 1 and tvA = 0). See |2C 
the evolution equations, and 
dicates why the BSSN formulation should be superior to 
ADM at least for linearized perturbations of flat space. 

In order to obtain the stable evolutions described be- 
low, we have found it necessary to add the following in- 
gredients to the BSSN formulation: 



1. As discussed in |29|, we actively force the trace of 



the conformal-traceless extrinsic curvature Aij to 
remain zero during our simulations by subtracting 
it after each time step. 

2. We use the independently evolved "conformal con- 
nection functions" F* only in terms where deriva- 
tives of these functions appear. Whenever these 
functions are undifferentiated, we recompute them 
from the conformal Christoffel symbols. We have 
found this to be very important to achieve long- 
term stability, but at the moment we lack a theo- 
retical understanding as to why this is so. 

Slicing conditions — As a first approach to evolving 
the solution described above, one could think of using 
the exact value of the lapse. It turns out that it is dif- 
ficult to keep the evolution stable if the lapse is not al- 
lowed to adapt to the (numerically induced) evolution 
of other dynamical quantities, particularly the trace of 
the extrinsic curvature. In order to obtain stable evo- 
lutions we have found it crucial to use a "live" slicing 
condition. What is required is a slicing condition that is 
well adapted to the exact solution in the sense that for 
this solution it recovers the exact lapse. For this we start 
from the Bona-Masso family of slicing conditions [pO| 



dta = —o? / [a) trif , 



(2) 



with /(a) > 0. As it is, this condition does not reproduce 
our exact solution for which \xK ^ 0, but dtOi=^. How- 
ever, one can easily see that for zero shift Eq. (||) implies 
dtOi cx 9t(det(7). For this to hold also with non-zero shift 
Eq. (0) must be generalized to 



dta = -af {a) [a tiK - \/i(3'] 



(3) 



For any static solution Eq. (^ implies dta=0. 

Another natural slicing condition to consider is 
dttiK=0. For initial data with tiK=0 this condition 
leads to maximal slicing, but dttrK=0 is a gauge choice 
that can be made in general, even if trK ^ 0, as is the 
case for the constant time slices of the black hole in EF 
coordinates. This "K freezing" condition leads to an el- 
liptic equation for the lapse, 



(4) 



In the numerical implementation, we solve this equation 
for the lapse but we hold tiK constant in time by hand. 
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In |g7|] in the context of the evolution of strong waves we 
have found that otherwise a drift away from the initial 
value due to numerical errors can lead to an instability. 
Such drifts were one of the reasons that led us to consider 
trace-split formulations like BSSN, because here trK is 
evolved as an independent variable which makes it trivial 
to enforce dttrK = 0. 

Shift conditions — In contrast to the experience with 
the lapse, we have found that using a static (exact) shift 
does allow us to get long-term stable evolutions. How- 
ever, this is not useful in general, so we have considered 
also live shift conditions. Live shifts have been studied 
before for black hole spacetimes in where a mini- 
mal distortion shift condition led to limited stability 
{t lOOM) for a single excised black hole. 

In our case a good choice was a conformal version of the 
3-harmonic shift ||3^ . 3- harmonic shifts play a natural 
role in mixed elliptic- hyperbolic systems p3t . The condi- 
tion we impose in the BSSN system is dtr'^=0 ("Gamma 
freezing" condition, note that T'' ^0), or 

(5) 

As mentioned before, 9fcr* is computed from the inde- 
pendent variable F*, in other terms we use F* = f ^'F^j, 
(notice that the momentum constraint was used to re- 
place djA'^ in d)). Equation (|) is an elliptic equation 
for the shift vector. For the solution of (^ and (||) we 
have used the multi-grid solver from BAM As in 

the case of K freezing, we explicitly hold the value of F* 
constant in time in order to prevent this quantity from 
drifting due to numerical errors. As shown in Section IV, 
allowing F* to drift results in an unstable evolution. 

We have also looked at shift prescriptions given by 
evolution equations instead of elliptic conditions. One 
way to do this is to transform an elliptic equation into a 
parabolic one by making 9t/3* proportional to the given 
elliptic operator ("driver" conditions, see Q). As an ex- 
ample we considered the following evolution equation for 
the shift obtained from the Gamma freezing condition (a 
"Gamma driver" condition) 

dtP'^kdtT', (fc>0). (6) 

Boundary conditions — There are two very different 
boundaries to consider in our simulations: the outer 
boundary of the numerical grid, and the inner bound- 
ary of the excised region. In principle there should be a 
rigorous treatment of numerical boundaries at finite radii 
(starting e.g. from js^, the first analytic treatment of the 
initial boundary value problem) . Here we are looking for 
simple numerical methods that are sufficient for the evo- 
lution of excised black holes. 



At the outer boundary we have attempted to keep all 
fields equal to their exact values, but have found that this 
introduces late time instabilities. Using a live boundary 
condition allows us to eliminate these instabilities. The 
boundary condition we use is a radiative boundary con- 
dition applied to the difference between a given variable 
and its exact value: / — /exact = u{r — t)/r. We apply 
this condition to all fields (even to the lapse and shift in 
the case of the algebraic gauge conditions) except the F' 
which we leave fixed to their exact values at the bound- 
ary. Applying this condition to the F* causes a drift away 
from the exact solution that eventually crashes the simu- 
lation (the origin of this drift is not well understood, but 
it seems to be related to the shift choice and is not present 
if one uses the Gamma driver shift described above). 

As to what boundary condition to use at the sides of 
the excision cube, we have experimented with many dif- 
ferent conditions and have finally settled on one that sim- 
ply copies the time derivative of every field at the bound- 
ary from its value one grid-point out along the normal 
direction to the cube (at edges and corners we define the 
normal direction as the diagonal). This condition is per- 
fectly consistent with evolving a static solution, where 
the time derivatives are supposed to be zero. Even in 
a dynamical situation, this condition is still consistent 
with the evolution equations since it is equivalent to just 
calculating the source term one grid point away. This 
means that our boundary condition should introduce a 
first order error, but as mentioned above, we do not ex- 
pect this error to affect the solution outside the horizon. 
One could in principle argue that nothing prevents gauge 
modes and constraint violating modes from propagating 
outside the horizon, thus spoiling the second order con- 
vergence of the exterior scheme. We have looked carefully 
at the convergence of our simulations, and have found no 
evidence that this happens in practice. 

IV. NUMERICAL RESULTS 

We now present some results of our numerical simu- 
lations. As discussed above, our simulations have been 
done with a live lapse condition, and we have considered 
both a static shift, and several live shift conditions. In 
our runs we have always taken M=l, so the horizon is 
a sphere of radius r=2, and we excise a cube of side 1 
(we have in fact also excised cubes of different size, but 
the results discussed below are not affected by this). The 
numerical integration is carried out using the so-called it- 
erative Crank-Nicholson scheme with 3 iterations (count- 
ing the initial Euler step as iteration 1). Because of the 
spherical symmetry of the problem typically only one oc- 
tant was evolved (with positive x, y, and z). However, 
as discussed at the end of this section, an unstable mode 
appears when the same simulations are performed on the 
corresponding full grids. 

Static shift — We first consider the case when the shift 
remains equal to its exact value. Figure l] shows a log 
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FIG. 1. Log plot of r.m.s. of the change in the lapse; 
Aa; = 0.4, At = 0.1. a) 53^ grid points, boundary at 20M. 
b) 103^ grid points, boundary at 40M. 



plot of the root mean square (r.m.s.) of the change in 
the lapse between consecutive time steps for two simula- 
tions using slicing condition with / = 1/a ("l+log" 
slicing [^,^), a grid spacing Aa; = 0.4, and a time 
step At = 0.1. Figure |l|a shows the results of a simula- 
tion using 53^ grid points, with the outer boundaries at 
20M. The change in the lapse drops as an exponentially 
damped oscillation until at t 3500M it reaches the 
level of round-off error (10~^^) and settles down (other 
functions show a similar behavior). The evolution was 
stopped a,t t — 4000Af, but it is clear that it could 
have continued. Figure |^b shows a simulation with the 
same resolution, but using 103^ grid points, with the 
outer boundaries now at 40M. The simulation goes past 
t ~ 2000M, and seems to have settled on an exponen- 
tially decaying oscillating pattern. (This simulation took 
100 hours running on 16 processors of an Origin 2000 
SGI machine. If the pattern continues, round-off error 
level would be reached by t ~ 12000M, requiring an- 
other 500 X 16 hours of computer time). The most ob- 
vious differences between the run with the boundaries at 
20Af and that with the boundaries at 40M is the fact 
that the period of the oscillations increases and the rate 
of decay decreases. The period increases by a factor of 
3.4 as we double the distance to the outer boundaries, 
so the oscillation time scale is not given directly by the 
light travel time from the boundary (which would ap- 
proximately double) . We do not know exactly what fixes 
this time scale, but the fact that when we look at in- 
dividual metric components we see that the oscillations 
behave like standing waves (and not travelling pulses) 
would seem to indicate that we are looking at different 
modes of oscillation of the whole system (interior plus 
boundaries). 

These simulations are not only stable for very long 
times, they are also exceedingly accurate. We have lo- 
cated the apparent horizon every 50 time-steps (using 
the 3D finder described in [p8[ ) , measured its area A and 
computed its mass M = (167r). Figure ^ shows the 
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FIG. 2. Evolution of horizon mass for the same simulations. 
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FIG. 3. Late time Hamiltonian constraint for runs with 
different resolutions. The values for the higher resolution runs 
were multiplied by factors of 4 and 16. 

behavior of the horizon mass as a function of time. In 
both cases, after an initial transient, the mass settles on 
a stationary value with an error of less than 1%. 

In Figure ^ we consider the convergence of our simula- 
tions by looking at the late time value of the Hamiltonian 
constraint along the x axis for simulations with 28^, 53'^ 
and 103"^ grid points and resolutions of Aj:=0.4,0.2,0.1 
respectively (boundaries at IQM). The Hamiltonian con- 
straint for the higher resolution runs has been multiplied 
by factors of 4 and 16. The fact that the three lines 
coincide indicates second order convergence. 

Elliptic shifts — We now consider results with elliptic 
shifts, such as those that we expect will be needed in a 
3D black hole merger simulation. Figure ^ shows two 
stable and three unstable runs up to t = 400M, and 
Figure || shows those three runs that lasted longer up to 
t = 3000Af. Second order convergence has been checked 
using two grids with 19^ and 35'^ points with the outer 
boundary at 7AI . For l-flog slicing a radiative boundary 
condition is applied to the lapse, while lapse and shift for 
the elliptic conditions are held fixed at the exact values. 

Stable runs are obtained for Gamma freezing shift with 
either l-flog or K freezing slicings. Referring to Figures 
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time (in M) 

FIG. 4. Log plot of r.m.s. of the change in the lapse for 
different lapse and shift combinations involving elliptic con- 
ditions; Aa; — 0.4, At = 0.1, 35"^ points, boundary at 13M. 
Run 1: stable (F freezing without drift, l+log); run 2: stable 
(F freezing without drift, K freezing without drift); run 3: 
unstable (F freezing without drift, l+log, static outer bound- 
aries); run 4: unstable (F freezing with drift, l-flog); run 5: 
unstable (minimal distortion, l-flog). 



§ and [HI for l+log slicing Aarms falls below 10^^^ at 
t ~ 1500M after four oscillations (run 1), while for K 
freezing there are more than fifteen oscillations, which 
damp out at around 10~^° followed by a straight line 
decay (run 2). 

The l+log. Gamma freezing run becomes unstable if 
the boundary values of all fields are static (run 3, crashing 
at t 1500M, Figure |b), or if dtt^ = is not set to zero 
identically and is allowed to drift because of numerical 
errors (run 4, crashing at 375M). We also tested l+log 
slicing with a minimal distortion shift | pT[ | computed from 
the ADM variables, but this run fails already at 27 M (run 
5). 

Algebraic shifts — Finally, we consider a simulation us- 
ing l+log slicing and a Gamma driver shift with k = 0.1. 
Figure |6| shows the r.m.s. of the change in the lapse 
and the horizon mass for a simulation with Ax = 0.4, 
At = 0.1 and 53^ grid points. After t ~ 2500M the 
solution becomes static up to round-off error. 

Discussion — The above results demonstrate that sta- 
ble 3D black hole runs can be obtained with the sim- 
ple excision technique that we introduced in this paper, 
with a variety of different gauge conditions. However, 
repeating these runs on a full grid as opposed to just 
one octant, with otherwise identical parameters, uncov- 
ers an unstable mode. Figure |^ shows as an example the 
situation for l+log slicing and static shift, although the 
problem appears for all the gauge conditions considered 
here. Tracing the growth of the unstable mode back in 
time suggests that it has started as numerical round-off 
error of around lO^^"' at t = 0. Increasing the grid reso- 
lution appears to have no significant effect on the growth 
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FIG. 5. Runs 1, 2, and 3 of Figure ^for run times of up to 
t = 3000M. 
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FIG. 6. Simulation using Gamma driver shift with fc=0.1; 
Aa; = 0.4, At = 0.1, 53^ points, boundary at 20A/. a) Log 
plot of r.m.s. of change in the lapse, b) Horizon mass. 



rate of the unstable mode, but the simulation now crashes 
slightly sooner. However, we do see good second order 
convergence at early times, before the instability becomes 
apparent. The situation does not improve if we impose 
the exact data at the excision boundary (imposing exact 
data at the excision boundary in octant mode works well 
and leads to stable simulations). Also, the presence of 
a horizon does not seem to be the cause of the problem 
since when we excise a cube that contains the horizon, 
as opposed to being contained by it, the instability is 
still present although it becomes somewhat milder (not 
surprising since we have excised a region with stronger 
data). While the achievable run times of about 500M 
are roughly 10 times larger than for singularity avoid- 
ing slicings, we have found that introducing an artificial 
asymmetry on the full grid by simply off-setting the ex- 
cision box one grid point in all directions makes the runs 
fail much sooner. Although the slope of the blow-up is 
not significantly affected when this artificial asymmetry 
is introduced, the exponential growth becomes evident 
from the very beginning. On the other hand, the full 
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FIG. 7. Unstable mode on a full grid for 1+log slicing with 
a static shift. Shown is a log plot of the r.m.s. of the change 
in the lapse for an octant run with Aa; = 0.4, At — 0.1, and 
28^ grid points together with the corresponding full grid run. 



grid runs can be stabilized by setting certain terms in 
the BSSN equations to their analytic values. In particu- 
lar, freezing the evolution of the F* while keeping the shift 
static suffices to obtain stability. In conclusion, the insta- 
bility appears to be more directly linked to the system of 
evolution equations than to the boundary condition, and 
we will investigate different variations of the evolution 
system in the future. 

We have also repeated the above simulations using the 
ADM equations with the same gauge and boundary con- 
ditions, and the same numerical techniques, but these 
runs fail already at i ~ 30M even in octant mode. 

V. CONCLUSIONS 

We have described a black hole excision technique in 
3-1-1 numerical relativity that has allowed us to obtain ac- 
curate, long-term stable evolutions of a black hole space- 
time in 3D. The main limitation is that the transition 
from octant symmetry to full grids introduces an unsta- 
ble mode, which is currently under investigation. Our 
implementation of excision is based on the idea of simpli- 
fying all ingredients of the excision algorithm as much as 
possible. In our case this means: 1) excising a cube nat- 
urally adapted to the underlying Cartesian coordinates, 
2) imposing a simple but stable boundary condition on 
the sides of this cube, and 3) using an upwind scheme 
instead of causal differencing. Crucial for obtaining our 
long-term stable evolutions has been the use of a live slic- 
ing condition and a radiation outer boundary condition. 
Although keeping a static shift does not appear to have a 
detrimental effect on the stability of our simulations, we 
have also experimented with several live shift conditions, 
both algebraic and elliptic, that can be generalized to 
more interesting physical situations. We consider these 
results a necessary first step towards the development of 
excision techniques capable of evolving the full inspiral 
collision of two black holes in an accurate and stable way. 
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